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Radiation driven outflow in active galactic nuclei: the 
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ABSTRACT 

We perform time-dependent, two-dimensional, hydrodynamical, numerical simu- 
lations to study the dynamics of a slowly rotating accretion flow from sub-pc to pc 
scales under the irradiation from the central AGN. Compared to previous work, we 
improve the calculation of the radiative force due to X-rays. More importantly, in 
addition to radiative pressure and radiative heating/cooling directly from the central 
AGN, in the momentum equation we also include the force due to the scattered 
and locally produced photons. We find that the accretion flow properties change 
significantly due to this "re-radiation" effect. The inflow rate at the inner boundary is 
reduced, while the outflow rate at the outer boundary is enhanced by about one order 
of magnitude. This effect is more significant when the density at the outer boundary 
is higher. The properties of outflows such as velocity, momentum and energy fluxes, 
and the ratio of outflow rate and the accretion rate, are calculated. We find that 
the efficiency of transferring the radiation power into the kinetic power of outflow is 
typically fO""^, far below the value of ^ 0.05 which is assumed in some cosmological 
simulations. The effect of the temperature of the gas at the outer boundary (Tq) is 
investigated. When Tq is high, the emitted luminosity of the accretion flow oscillates 
because of the strong radiative heating. Another question we hope to address by this 
work is the so-called "sub-Eddington" puzzle. That is, observations show that the 
luminosity of almost all AGNs are sub-Eddington, while theoretically the luminosity 
of an accretion flow can easily be super-Eddington. We find that even when the 
re-radiation effect is included and outflow does become much stronger, the luminosity, 
while reduced, can still be super-Eddington. Other observational implications and 
some caveats of our calculations are discussed. 

Key words: accretion, accretion discs - hydrodynamics - methods: numerical - galax- 
ies: active - galaxies: nuclei 



1 INTRODUCTION 

Active galactic nuclei (AGNs) are thought to play an im- 
portant role in the processes of galaxy formation and evolu- 
tion. The most important evidence is perhaps the remark- 
able correlations between host galaxy properties and the 
mass of the supermassive black holes (e.g., Magorrian et 
al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; 
Graham et al. 2001; Kormendy et al. 2009) observed in the 
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past decade. Many theoretical studies have also confirmed 
that AGN feedback can affect the ambient environment ef- 
fectively (e.g.. Silk & Rees 1998; Ciotti & Ostriker 1997, 
2001, 2007; Ciotti, Ostriker & Proga 2009; King 2003; Mer- 
loni et al. 2004; Springel, Di Matteo & Hernquist 2005a; 
Sazonov et al. 2005; Di Matteo, Springel, & Hernquist 2005; 
Hopkins et al. 2005; Murray, Quataert & Tompson 2005; 
Somerville et al. 2008; Ostriker et al. 2010). 

The AGN feedback can be in the forms of radiative (e.g., 
Ciotti & Ostriker 1997, 2001, 2007, Ciotti et al. 2009) or me- 
chanical energy input (e.g., Springel, Di Matteo & Hernquist 
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2005b; Ostriker et al. 2010). In terms of the radiative feed- 
back, both momentum feedback and energy feedback work, 
i.e., the ambient gas could be accelerated by radiation pres- 
sure and heated by radiative heating. Mechanical feedback is 
via wind/outflow. For example, in the context of hot accre- 
tion flow (e.g., advection-dominated accretion flow). Yuan, 
Xie & Ostriker (2009; see also references therein) studied 
the effect of global Compton scattering, namely the scatter- 
ing between the photons produced at the innermost region 
of the ADAF and the gas at large radius. They found that 
when the luminosity is higher than ~ 2%LBdd, the heating 
is so strong that the flow will be heated to be above the 
virial temperature thus the accretion will be stopped. It is 
shown that this physical mechanism is able to explain the 
intermittent activity of some compact young radio sources 
(Yuan & Li 2011). 

Provided the gas is moderately ionized and can interact 
with the UV continuum through many UV line transitions, 
radiation pressure force due to spectral lines (i.e., line force) 
has been shown to be very efficient at producing winds from 
accretion disks. This can occur in the small scale of the 
innermost accretion flow (e.g., Murray et al. 1995; Proga, 
Stone & Kallman 2000, hereafter PSKOO) and much larger 
scale from ~ 0.01 pc to 10 pc (Proga 2007a; Proga, Os- 
triker & Kurosawa 2008; Kurosawa & Proga 2009, hereafter 
KP09). Some highly blueshifted broad absorption line fea- 
tures seen in the UV and optical spectra of AGN can be 
explained by this mechanism (e.g. Chartas, Brandt &: Gal- 
lagher 2003; Crenshaw, Kraemer & George 2003; Proga &: 
Kallman 2004; Hamann et al. 2008). However, not all AGN 
outflows can be explained by line-driven mechanism because 
of, e.g., the too-low luminosity, too-high ionization state, or 
both (e.g., Chelouche & Netzer 2005; Kraemer et al. 2006; 
see Proga 2007b for a review). This has also been shown to 
be the case of black hole X-ray binaries which are believed 
to be the scaled-down version of AGNs (e.g.. Miller et al. 
2006, 2008; Luketic et al. 2010). 

Therefore, other mechanisms are required, such as ther- 
mal (e.g., Begelman, de Kool & Sikora 1991) and, espe- 
cially, magnetic driving (e.g., Blandford & Payne 1982; Em- 
mering, Blandford & Shlosman 1992). Most recently, based 
on global HD and MHD numerical simulations. Yuan, Bu 
& Wu (2012; see also Li, Ostriker & Sunyaev 2013) show 
that outflow must exist in hot accretion flows. Moreover, 
Yuan et al. proposed a Blandford & Payne-like mechanism 
of producing winds. This mechanism should be quite uni- 
versal since it is a result of magneto-rotational instability 
(MRI), which is widely believed to exist in accretion flows 
and is the mechanism of transporting the angular momen- 
tum. Previous works on studying the production of MHD 
winds usually require the pre-existence of large scale poloidal 
fleld and treat the disk as the boundary condition. Different 
from these works. Yuan et al. (2012) simulate the disk and 
outflow self-consistently and simultaneously, and the pre- 
existence of the large-scale poloidal magnetic field is not re- 
quired. However, the simulation was done for hot accretion 
fiows. Thus the next step it remains to be probed explicitly 
whether it can be equally effective for thin disks, although 
in principle it is expected to work. 

Among the above-mentioned line-driven works, KP09 
studied the momentum and energy interaction between the 
radiation coming from the central AGN and the accretion 



fiow in the region between ~ 10^^ pc and 7 pc. They focus 
on the line-driven outfiow. Compared with previous works, 
the luminosity of the central AGN is not fixed, but self- 
consistently determined by evaluating the mass accretion 
rate at the inner boundary. They found that outflow can be 
driven efficiently from the infiow. When the temperature of 
the gas at the outer boundary is high. To > 2x lO'^K, a strong 
correlation between the mass outflow rate (Mout) at the 
outer boundary and the luminosity (L) exists, Mout oc L'^ 
with g ~ 2. When To = 2 x lO^K and when L > LEdd, 
the correlation becomes flatter, < g < 0.2. Another in- 
teresting result is, they found that when the density at the 
outer boundary is high, the accretion luminosity of the cen- 
tral AGN can well exceed the Eddington value. The super- 
Eddington accretion mainly proceeds in the equatorial plane 
and is possible because the radiation flux from the disc is 
signiflcantly reduced in the equatorial direction due to the 
geometrical foreshortening effect. 

In almost all previous works on radiatively driven out- 
flow including KP09, they properly consider the forces due 
to electron scattering and UV line processes due to the "first- 
order" radiation coming from the central AGNs. However, 
the effect of the locally produced photons via scattering, 
the bremsstrahlung, and line radiation, are all neglected. 
Thus when radiation is absorbed and then re-emitted the 
dynamical effects of the latter are ignored, essentially vio- 
lating energy conservation. In principle these photons could 
play an important role. For example, they can produce an 
additional force in the vertical direction and thus puff the 
disk up. This may make the interaction between the radia- 
tion and the gas significantly stronger. In the present paper 
we want to investigate this "re-radiation" effect. 

Another motivation to consider this effect in the present 
work is that we hope to address the so-called following 
"sub-Eddington" problem. Observations to a large sample 
of AGNs with 407 sources show that the distribution of their 
Eddington ratio (= Lboi/^Edd) is well described by lognor- 
mal, with a peak at Lboi/iEdd ~ 1/4 and a small dispersion 
of 0.3 dex (KoUmeier et al. 2006). In other words, almost all 
AGNs are radiating below LEdd- Later Steinhardt & Elvis 
(2010) extended this study to a much larger sample con- 
sisting of 62185 quasars from the Sloan Digital Sky Survey 
and they reached the similar conclusion (but see Kelly & 
Shen for a different opinion). On the other hand, it has long 
been known from the black hole accretion theory that the 
rotating accretion fiow can well radiate above Lsdd- Such an 
accretion fiow is called slim disk or optically thick advection- 
dominated accretion flow (Abramowicz et al. 1988). This 
analytical result has later been confirmed by radiative hy- 
drodynamical numerical simulations (Ohsuga et al. 2005). 

It has been argued that the galactic-scale fueling events 
are likely to have a broad mass distribution so this should 
not be the reason for the sub-Eddington luminosity of AGNs. 
One possible way to solve this puzzle is that, the accre- 
tion rates are determined by the self-regulation of the ra- 
diative feedback. In other words, the strong radiation from 
the inner accretion fiow may drive a strong outfiow thus 
reducing the accretion rate to below the Eddington rate 
{M-Edd = lOLEdd/c^). This should not be effective in the 
innermost region close to the black hole, since the interac- 
tion between radiation and accretion flow has been properly 
considered in Ohsuga et al. (2005) . At larger scale, however. 
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the simulations by KP09 still find super-Eddington accretion 
even though the radiative feedback is included in their simu- 
lations. In this paper, we include the re-radiation effect to do 
the calculation again, to see whether the "sub-Eddington" 
puzzle can be solved. It is vitally important that the effects 
of radiation feedback at the Bondi radius be correctly in- 
cluded, since it is in this region that the rate of inflow is 
determined. 

In this paper, following the approach of KP09, 
we perform two-dimensional hydrodynamical (HD) time- 
dependent simulations of accretion flow irradiated by the an- 
gle dependent UV flux from an accretion disk and isotropic 
X-ray flux from a corona. The models include radiation 
forces due to electron scattering and line processes from both 
the first-order radiation and the re-radiation, and radiative 
cooling/heating. Special attention is paid to the re-radiation 
effect. In §2, we describe some basic aspects of our model, es- 
pecially the radiative force due to re-radiation. The numer- 
ical simulation method and the initial and boundary con- 
ditions are also described in this section. The results of the 
simulations are given in §3, and §4 provides some discussion. 
Finally, §5 is devoted to a summary. 



2 MODEL 

The basic scenario of the model, together with most of the 
formula for calculating the radiative cooling and heating, are 
the same with those presented in KP09 (see also PSKOO). 
The readers are referred to that paper for details. However, 
we do make some changes compared to KP09. In addition to 
the inclusion of reradiation force mentioned in §1, we also 
improve, we believe, the calculation of the radiation force 
due to X-ray. 



2.1 Hydrodynamics of the accretion flow 

Our aim is to study the dynamics of the accretion flow be- 
tween an inner radius r^ and outer radius Vo- In addition to 
the gravitational force, this gas will inevitably be irradiated 
by the radiation emitted from the central AGN. The radia- 
tion will heat and push the gas by Compton scattering and 
line absorption, if the gas is not fully ionized. The dynam- 
ics of the accretion flow is described by the following set of 
equations. 



Dp 
Dt 



-f pV ■ V = 0, 



p^ = -vp+pg+pF-^ 



_D 



-pV ■ V + pC, 



(1) 



(2) 



(3) 



where p, e, p, and v are the mass density, internal energy 
density, gas pressure, and velocity, respectively; p£. describes 
the radiative cooling and heating, simply as the net heating 
rate (see PSKOO) ; g is the gravitational acceleration of the 
central object; F"^'^'^ (see equation I12p is the total radiation 
force per unit mass. We adopt an adiabatic equation of state 
p = (7 — l)e, and consider models with the adiabatic index 
7 — 5/3, with, of course, allowance for heating and cooling 



via the C term in equation (3). The implementations of F'^'''^ 
and pC are described as follows. 

2.2 The central AGN 

The central AGN is described by a super-massive black hole 
with mass Mbh surrounded by an accretion flow. The na- 
ture of the accretion flow in a luminous AGN is still an un- 
solved problem. Usually people assume that it is described 
by a standard thin disk. However, a thin disk alone can't 
explain the origin of hard X-ray emission widely detected in 
AGNs. For simplicity, we assume that the accretion flow is 
a combination of a standard thin disk (Shakura & Sunyaev 
1973) plus a spherical corona. They are responsible for the 
optical/UV and X-ray emissions respectively. The radiation 
from the corona is assumed to be isotropic. 

In the point-source approximation limit, the radiation 
flux from the central X-ray corona can be written as 

/*I/accexp(— rx) 



T, = 



47j-j,2 



(4) 



where /* is the fraction of the total luminosity (I/acc) in X- 
ray. We adopt /« — 0.05 throughout the paper, rx is the 
X-ray optical depth in the radial direction. 



TX 



pKx dr, 



(5) 



where kx is the X-ray opacity. We assume that the atten- 
uation is dominated by Thomson scattering. We usually 
set Kx = 0.4 cm g^ . However, this must over-estimate 
the scattering opacity. This is because unlike true absorp- 
tion, scattering merely re-directs the photons. The scattered 
X-ray photons must be replenished by photons scattered 
from other propagation lines. So another extreme is to take 
Kx = 0. In the present paper, we also test this case for 
comparison. The real opacity should be a certain value be- 
tween these limits. The radiation from the standard thin 
disk J-d ex |cosS|. Consequently, the disc radiation flux at a 
distance r from the center can be written as 



Td = 2 jcos^l 



fdLaccexp{—Tuv) 

47f}-2 



(6) 



where /d = 1 — /» is the fraction of the total luminosity in 
the disc emission, tuv is the UV optical depth and we as- 
sume Tuv ~ 0. Note that the "|cosS|" assumption is poten- 
tially a strong assumption, since it suppresses the radiation 
force close the equatorial plane if we include re-radiation. 
We therefore have done several test simulations by replac- 
ing "|COS0|" with "{Tperp + | COS 6*1 )/(2rperp + 1)"- Here Tprep 

is the scattering optical depth above the disk plane. We flnd 
that this effect does not influence our results significantly. 

In our simulation, at each time step, the accretion lumi- 
nosity will be calculated self-consistently based on the mass 
inflow rate at the inner boundary r^. It is given by 



Lace (t) = vM^{t)c^ 



(7) 



where Ma (t) is the mass inflow rate at r^ and at a given time 
t; rj is the radiative efficiency. We adopt rj = 1/12 in this pa- 
per. Following KP09, the mass accretion rate at a given time 
t is actually time-averaged over some time interval to reduce 
the fluctuation level. Obviously, a time lag r is expected be- 
tween the change of this accretion rate to the disk and the 
change of the luminosity from the central AGN Lace. Since 
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in the present paper the accretion flow we consider has small 
angular momentum, a small disk with a radius of Rd will be 
formed within the inner boundary of our simulation. In this 
case, the time lag is roughly equal to the accretion timescale 
at the outer boundary of this small disk. If the small disk 
can be described by a standard thin disk, the lag time r can 
be approximated as tacc ~ Rd/vr = Rd/{aCaH/Rd); if it is a 
slim disk, the time lag should be tacc ~ Rd/vr ~ Rd/{avk)- 
Here a is the viscosity parameter, and we set a = 0.03 in 
this paper. The black hole mass is assumed to be constant. 
The value of t is not important for our aim. 

One caveat here is that we assume the accretion rate 
onto the black hole is equal to the inflow rate at r^. This 
may not be true, however. Almost all global numerical sim- 
ulations to hot accretion flows show that the accretion rate 
should decrease with decreasing radius (see Yuan, Wu & Bu 
2012 for a review). The reason is identified to be because of 
continuous mass loss in outflow produced by an MHD pro- 
cess similar to the Blandford-Payne mechanism (Yuan, Bu 
& Wu 2012). This mechanism in principle also works in the 
case of a cold thin disk. Actually, the initial global MHD 
numerical sinmlations to thin disk by Ohsuga & Mineshige 
(2011) do show significant outflow, although the exact mech- 
anism remains to be probed. Since the outflow in the context 
of a thin disk is still an open question, in the present paper 
we simply assume that the accretion rate is constant within 
ri and thus our computed luminosity is an upper bound on 
the true luminosity for given outer boundary conditions. 



2.3 Heating/cooling and radiative force of the gas 

The interaction between the radiation and the gas and 
the radiative processes we consider includes Compton 
heating/cooling. X-ray photoionization and recombination, 
bremsstrahlung, and line cooling. The net cooling rate de- 
pends on the photoionization parameter ^ and the character- 
istic temperature (or radiation temperature) of the radiation 
from the central AGN (Tx)- The gas photoionization state 
is determined by the photoionization parameter 



K = = j-tLaccexp(—Tx)/nr , 



(8) 



where n — p/^nip is the number density of the local gas 
located at distance r from the central AGN, /x is the mean 
molecular weight and is fixed to be 1. Only X-ray radia- 
tion ionizes the gas. The sum of photoionization heating- 
recombination cooling rate (n^Gx) and Compton heat- 
ing/cooling rate (n^Gcomp) is 

Ex — n {Gx + Gcomp), 



(9) 



with 



Gx = 1.5 xlO"^'^'/*r"''"(l-r/rx) ergcm''s-\ (10) 

Gcon.p = 8.9 X lO""'' C (Tx - T) erg cm^* s^V (11) 

Here Tx ~ 8 x lO^K (Sazonov et al. 2004) is the "character- 
istic temperature" of the X-ray radiation and it is 4 times 
larger than the Compton temperature (refer to PSKOO and 
Proga 2007a). For a low-luminosity AGN, Tx will be much 
higher, ~ 3 x lO^K (Yuan, Xie & Ostriker 2009). 



2.4 Radiative force 

Compared to KP09, in our present work we make two 
changes when calculating radiative forces. The radiative 
forces considered in KP09 are due to Compton scattering 
and line processes. The former is the only available radia- 
tive force if the gas is fully ionized. Obviously, this force is 
important only when the luminosity from the central AGN 
approaches the Eddington value. However, the force can be 
significantly enhanced if the gas is moderately ionized, since 
the opacity is much higher than scattering opacity due to 
the bound-bound and bound-free transitions. So we should 
include the line force. For this force, following KP09, we as- 
sume that only optical/UV photons have contributions, but 
neglect the possible force line due to some metal line in the 
soft X-ray band. Regarding the sources of photons, as we 
have emphasized in §1, in addition to the radiation from the 
central AGN, the "local" photons from the local radiative 
processes and scattered photons originally from the central 
AGNs are also included. Both of them contribute to the ra- 
diative force. 



F" 



r ' ^ z 



(12) 



The inclusion of the force due to the "re-radiation" process 
^prad,re-j j^ ^j^^ g^.^^ change we make compared to KP09. It 
turns out that the inclusion of this effect does significantly 
change the results, as we will state in the following section. 

For X-ray photons, KP09 only consider the force due to 
Compton scattering. However, if the gas is not fully ionized 
as in our present case, absorption processes such as pho- 
toionization should also contribute to the force. This force 
should be characterized by Ex/c. The inclusion of this force 
is the second change we make compared to KP09. 

Based on the above considerations, the total radiative 
force due to the central AGN photons then can be written 
as: 



Fr'*''= (r,0) = :^ + — 



i-'ac 



47rr2 



[l-f-Mllcos^l/d-f. (13) 



Here Kes = 0.4 cm^ g~^ is the scattering coefficient for free 
electrons, and M is the force multiplier (Castor et al. 1975) 
- the numerical factor which parameterizes by how much 
spectral lines increase the scattering coefflcient (see PSKOO 
for details) . We have calculated the dependence of the multi- 
plier M on the temperature T and the ionization parameter 
^. We found that the multiplier M could be up to the order 
of 10^ ~ 10'', providing that the gas is moderately ionized. 
It drops rapidly when temperature T or the ionization pa- 
rameter f is larger than some critical value (eg., ^ > 10, or 
T > 10^). Specifically, it is safe to ignore the line force (< 1% 
of electron scattering force) when ^ > 10'' or T > 2 x 10^, 
since then the gas is almost fully ionized. 

Now we calculate the second term in equation (|12p . i.e, 
the re-radiation force. Exact treatment of this force requires 
the full radiative transfer calculation which is beyond the 
scope of the present paper. Because the accretion fiow is 
rotating (see below), the fiow has a disk-like shape. In this 
case, for simplicity we can adopt the assumption of plane- 
parallel approximation. The "re-radiation" photons would 
eventually escape from the accreting fiow, and produce a 
net force in the vertical direction. Consider a rectangular 
box located within the disk, with the bottom of the box 
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Table 1. Summary of models with different parameters. 



Model 


Model 


ri 




PO 




To 


Re-radiation 


i^acc 


M 


,ut(r 


,) 


Vw 


Number 




(r.) 


(10- 


-21g |,jjj- 


-'■') 


(106K) 




(LEdd) 


(10= 


5g g- 


-') 




(1) 


(2) 


(3) 




(4) 




(5) 


(6) 


(7) 




(8) 




(9) 



1 


K5a 


500 


10 


2 


no 


1.19(0.19) 


15.31(9.12) 


0.76(9.12) 


2 


K6a 


500 


20 


2 


no 


1.85(0.43) 


25.61(24.29) 


0.84(0.92) 


3 


K6c 


1250 


20 


2 


no 


1.72(0.14) 


29.05(17.12) 


0.98(0.57) 


4 


K6d 


2550 


20 


2 


no 


1.71(0.23) 


31.17(13.22) 


1.10(0.53) 


5 


K7a 


500 


50 


2 


no 


3.97(1.10) 


21.54(27.46) 


0.35(0.58) 


6 


K7b 


625 


50 


2 


no 


3.69(0.78) 


20.80(17.05) 


0.33(0.24) 


7 


K7d 


2500 


50 


2 


no 


3.14(1.03) 


49.82(42.21) 


1.03(0.95) 


8 


K8a 


500 


100 


2 


no 


5.74(0.70) 


21.52(14.41) 


0.22(0.13) 


9 


K8d 


2500 


100 


2 


no 


5.70(1.75) 


34.24(50.46) 


0.35(0.51) 


10 


R5a 


500 


10 


2 


yes 


0.77(0.04) 


37.95(8.71) 


2.88(0.70) 


11 


R5b 


625 


10 


2 


yes 


0.94(0.04) 


31.53(5.07) 


1.96(0.34) 


12 


R5c 


1250 


10 


2 


yes 


1.10(0.10) 


24.28(10.67) 


1.28(0.55) 


13 


R5d 


2500 


10 


2 


yes 


1.12(0.11) 


26.88(8.05) 


1.40(0.45) 


14 


R6a 


500 


20 


2 


yes 


1.19(0.16) 


67.30(10.90) 


3.35(0.78) 


15 


R6b 


625 


20 


2 


yes 


1.28(0.17) 


65.14(10.73) 


3.02(0.69) 


16 


R6c 


1250 


20 


2 


yes 


1.48(0.18) 


39.10(12.13) 


1.54(0.45) 


17 


R6d 


2500 


20 


2 


yes 


1.60(0.18) 


44.04(11.45) 


1.61(0.45) 


18 


R7a 


500 


50 


2 


yes 


2.02(0.02) 


139.13(5.03) 


3.99(0.14) 


19 


R7b 


625 


50 


2 


yes 


2.03(0.03) 


137.86(4.83) 


3.93(0.14) 


20 


R7c 


1250 


50 


2 


yes 


1.81(0.03) 


127.04(11.80) 


4.08(0.39) 


21 


R7d 


2500 


50 


2 


yes 


1.60(0.03) 


194.29(1.92) 


7.04(0.14) 


22 


R8b 


625 


100 


2 


yes 


3.14(0.27) 


202.05(28.40) 


3.75(0.59) 


23 


R8c 


1250 


100 


2 


yes 


4.56(1.09) 


36.45(54.69) 


0.49(0.82) 


24 


R8d 


2500 


100 


2 


yes 


4.14(1.20) 


126.12(97.44) 


2.14(1.91) 


25 


R7C-X0 


1250 


50 


2 


yes 


1.76(0.22) 


120.61(9.67) 


4.03(0.58) 


26 


R7C-1T 


1250 


50 


0.5 


yes 


1.61(0.27) 


28.59(11.05) 


1.08(0.50) 


27 


R7C-1TX0 


1250 


50 


0.5 


yes 


1.19(0.01) 


52.17(3.28) 


2.54(0.16) 


28 


R7c-hT 


1250 


50 


6 


yes 


2.94(0.04) 


245.72(10.33) 


4.85(0.22) 


29 


R7c-hTX0 


1250 


50 


6 


yes 


2.26(0.46) 


338.04(29.36) 


8.98(1.85) 


30 


R7c-hT2 


1250 


50 


10 


yes 


2.02(0.97) 


461.93(100.00) 


15.43(5.07) 


31 


R7c-hT2X0 


1250 


50 


10 


yes 


2.74(0.51) 


508.50(57.05) 


11.18(3.12) 



Note: Values in brackets in Columns 7, 8 and 9 are the normalized standard deviations cr„ of the time series values. The symbol 'XO' 
in the model names means kx = in those models, 'IT' means lower temperature To, 'hT' means higher temperature Tq. 



overlapped with the equatorial plane of the accretion disk 
while the upper side at height z. "Re-radiation" photons are 
steadily produced within this box and then escape. Given 
the plane-parallel approximation and symmetry, the net ra- 
diative flux due to re-radiation (and thus the force) should 
be only in the upper side of the box. Applying the Gauss 
theorem, we can easily calculate the vertical radiative force. 



Fr ■'■"(r,e): 
where 



C Jo 



{Sc + n Ltrem + n Lli„e)dz ■ z, (14) 



-'-'ace 



[/*exp(-rx) + 2|cose|/d] 



(15) 



47J-J.2 

is the source term due to the first-order scattered photons 
of the radiation from the central AGN; while 



Lhicm = 3.3 X 10 T ' erg cm s 



(16) 



and 



Liine = s\i.7 X io"^*cxp(-i.3 X ioVr)r^r"^/^ 

(17) 
+ 10 Jcrg cm s 

are the rate of bremsstrahlung and line cooling, respectively. 



The parameter 5 is introduced to control line cooling, S = 1 
represents optically thin cooling, and S < 1 represents op- 
tically thick cooling. In general, the line cooling dominates 
over the other cooling processes. The formulae above only 
show the re-radiation force caused by electron scattering 
while the line force is ignored. So the effect of "re-radiation" 
calculated in this paper should be regarded as a lower limit 
in this sense. We did include the line force in some models. 
We found that the change of the result is not very significant. 
This is perhaps because the re-radiation force monotonically 
increases with height, and at a large height, the gas temper- 
ature becomes high so the line force is not important. We 
should emphasize that all our treatment is approximate, and 
a full radiative transfer calculation is required in the future. 



2.5 Simulation setup 

We solve the dynamical equations (1-3) using the code 
ZEUS-MP (Stone & Norman 1992a; Hayes et al. 2006). 
The simulations are performed in spherical polar coordinates 
(r, 6, <j}) assuming axial symmetry about the disc rotational 
axis (^ = 0° ) and in two-dimensions. Our computational do- 
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main is defined to occupy the angular range 0° ^ ^ 90° 
and the radial range ri ^ r ^ Tq. In KP09, both n and ro are 
fixed, i.e, they are same in all models. In the present work, we 
want to check the convergence for various values of ri , so we 
set various ri in different models but adopt ro = 2.5 x 10"'r,, 
where r* = Gvg is the innermost stable circular orbit (ISCO 
hereafter) of a Schwarzschild BH and rg = GM-bh/c? ■ We 
set the BH mass, Mbh ~ 10* M©. Our standard numerical 
resolution in the r direction consists of 144 zones with the 
zone size ratio, dvi+i/dri — 1.04 and consists of 64 zones 
in the 6 direction with d6j+i/d9j = 1.0. Gridding in this 
manner ensures good spatial resolution close to the inner 
boundary. 

For the initial condition, following KP09, we set the 
density and temperature of gas uniformly , i.e., p = po and 
T = To everywhere in the computational domain. The ini- 
tial velocity of the gas is assigned to have the distribution 
described by the following latitude-dependent angular mo- 
mentum at the outer boundary. 



m^ioii 



lo = VGMbht 



(18) 



where rdr is the "circularization radius" on the equatorial 
plane. In the present work, we set r^ir = 300r* < Vi. The 
boundary conditions are specified in the following way. We 
apply an axis-of-symmetry boundary condition at the pole 
(i.e., 9 = 0°) and reflecting boundary conditions at 6 = 90°. 
For the inner and outer radial boundaries, we apply an out- 
flow boundary condition (i.e., to extrapolate the flow beyond 
the boundary, we set values of variables in the ghost zones 
equal to the values in the corresponding active zones, see 
Stone & Norman (1992a) for more details). To represent 
steady conditions at the outer radial boundary, during the 
evolution of each model we apply the constraints that all 
hydrodynamical variables except radial velocity Vr (always 
let it float) in the last zone in the radial direction to be equal 
to the initially chosen values, i.e., p — po, T = To, ve = Q 
and v^ = I /{r sin 9), when Vr{ro,6) < (inflowing). We al- 
low all hydrodynamical variables to float when Vr{ro, 6) > 
(outflowing). This approach is to mimic the situation where 
there is always gas available for accretion. 

So the free parameters in our model include the inner 
radius Vi of the computational domain, the density po and 
temperature To of inflowing matter at the outer boundary. In 
this paper we focus mainly on the case of To = 2 x lO^K. But 
we do also consider other To for comparison. Various Vi and 
Po are adopted. The details of various models are listed in 
Table[T] The "yes/no" shown in Column 6 indicates whether 
re-radiation is included in the simulations. 



3 RESULTS 

Some results of various models are given in Table [l] Es- 
pecially, Column 7 gives the time-averaged luminosity cal- 
culated according to the mass accretion rate at the inner 
boundary. This may exceed the luminosity as seen at infln- 
ity due to absorption of radiation within the flow. Models 
1, 2, 5 and 8 are identical to the models 5, 6, 7 and 8 in 
KP09 except that we use our improved term when calcu- 
lating X-ray force, i.e., the flrst term in equation (|13|) . The 
last characters in the model name, i.e., "a", "b", "c", and 
"d", correspond to n = 500, 625, 1250, 2500 in units of r„ 



respectively. We deflne the mass inflow and outflow rates as 
below: 

Minir) = 27rr^ / p min(i;^,0) sin(6l) d9, (19) 

Jo 



Mout(r-) = 27rr^ / p max(i;^,0) sin(e) d9. (20) 

The time-averaged mass outflow rate is given in Column 8. 
The Column 9 gives 



Vw 



Ma ' 



(21) 



which is the ratio of mass outflow rate at the outer bound- 
ary and mass accretion rate at the inner boundary. It is 
important to note that the fraction of the material which is 
inflowing at the outer boundary which actually accretes on 
the black hole, deflned as fac = Ma/Min{ro), is from the 
mass conservation fac = 1/(1 + 7?m)i ^o that in cases where 
t;™ 2> 1 the net accretion, face is small. 



3.1 Effects of re-radiation 

We take models 7 series (i.e., K7* and R7*) as our fiducial 
models and explore their properties in details. Fig. [T] com- 
pares the radial profiles of mass infiow and outfiow rates, 
gas density and temperature between models K7a and R7a. 
The former does not consider re-radiation while the latter 
does. Fig. [2] shows the snapshots of contours of gas density 
and temperature of the two models. From the two figures 
we can find the following changes after re-radiation is taken 
into account. 

(i) Close to the equatorial plane the gas density decreases 
when re-radiation is included, as shown by the bottom-left 
panel of Fig. [T] However, we can see from Fig. [5] that away 
from the equatorial plane the density becomes higher, i.e., 
the density scale-height becomes larger. The reason is that 
the accretion fiow expands upward due to the inclusion of the 
vertical radiative force from re-radiation photons, then the 
vertically expanded inner edge of the accretion fiow shields 
the radiation from the central AGN, which makes the den- 
sity higher and the mass infiow rate larger. 

(ii) From Fig.[2l we can see that the temperature becomes 
lower in most of the region. This is because the increase of 
density produces stronger radiative cooling. The decrease 
of temperature close to the equatorial plane is because the 
compression work induced by accretion becomes weaker due 
to the inclusion of the vertical force. 

(iii) Perhaps the most important effect of including re- 
radiation is that outfiow becomes nearly one order of mag- 
nitude stronger, as clearly shown by the comparison between 
the top two panels of Fig. [l] This is easy to understand. As 
we state above, the inclusion of the vertical force due to 
re-radiation photons makes the density scale-height larger. 
Since the irradiation force from the central AGN oc cos 9, 
the gas will be more easily blown away if located at higher 
latitude. 

(iv) The mass inflow rate at the inner boundary decreases 
by a factor of 2. This is obviously because of the enhanced 
mass loss via the outflow. It is interesting to note that the 
mass inflow rate at the outer boundary also increases by a 
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Figure 1. Comparison between models Model K7a (without considering re-radiation; top left panel and the red solid line in the bottom 
panels) and Model R7a (with re-radiation included; top right panel and the blue dotted line in the bottom panels). Top panel — The 
radial profiles of mass inflow rate (blue dotted line), outflow rate (green dashed line), and the net rate (red solid line) in unit of lO'^^g s~^ . 



Bottom panel — The radial profiles of gas density {left panel, in unit of 10 
over three grids above the equatorial plane. 



g cm ) and temperature (right panel, in unit of K) averaged 



factor of 2. The reasons are two-fold. One is that the radia- 
tion from the central AGN decreases because of the decrease 
of the inflow rate at the inner boundary, thus the outward 
radiation force becomes weaker. Another reason is that, as 
we stated above, the gas density in most of the region in- 
creases. This gas will more effectively shield the gas at the 
outer boundary from irradiation by the central AGN. 

(v) The flow morphology also changes with the inclusion 
of re-radiation. Comparing the top and bottom panels of 
Fig. O we can see that in Model K7a, the inflow falls inward 
in a relatively small range of polar angle, i.e., 9 ~ 80° — 90°. 
But in Model R7a, i.e, the re-radiation effect is included, 
the inflow occurs in a relatively larger range, 60° < S < 90° . 
In addition, it is interesting to note that above the inflow 
region, the motion of the gas is turbulent in Model K7a; 
while in Model R7a, the motion is less turbulent. 

KP09 studied the correlation between the mass out- 
flow rate at the outer boundary and the luminosity in 
units of Eddington luminosity r(= L/Z/Edd) where Z/Edd ~ 
AncGMBii/o'e) is the Eddington luminosity. They found that 
when To is not too high. To < W^K, outflows exist only 



when r > 0.5. Speciflcally, when To = 2 x l(fK, as in our 
paper, a strong correlation is found for 0.5 < F < 1.5, i.e., 
Mout oc F' with q = 2.3(±0.3). The correlation index q only 
becomes slightly smaller when the temperature To at the 
outer boundary is higher. In the case of To = 2 x 10^ K, Mout 
becomes almost saturated when F > 1.5. KP09 explains this 
result as because the inflow is squeezed to a narrow range 
of when F becomes high. 

We have studied the effect of re-radiation on the corre- 
lation. Fig.Oshows the result. All the lines and dots except 
the three solid blue dots are taken from KP09. The three 
solid blue dots are the results when re-radiation effect is in- 
cluded for To = 2 X lO^K. We can see that the correlation 
index q roughly remain the same, and the correlation again 
becomes almost saturated when F > 2. The only change is 
the normalization: the outflow rate increases by almost one 
order of magnitude as we have stated above. 

At last, as we have mentioned above, the re-radiative 
force makes the density in most of the region larger. In this 
case, the attenuation of X-rays also becomes much more 
significant, thus the photoionization parameter ^ decreases. 
As we have mentioned in section [2^ the line force multiplier 
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Figure 2. Snapshots of contours of logarithmic gas density (left column) and logarithmic temperature (right column) at steady state. 
The top panel is for Model K7a (without re-radiation), and the bottom panel is for Model R7a (with re-radiation included). 



usually increases with decreasing ^. In other words, the line 
force should increase rather than decrease when we consider 
the re-radiation effect. 



3.2 Convergence with the inner boundaries 

It is necessary to check whether the result converges with 
various inner boundaries r^. With this aim, we have imple- 
mented a series of runs (see Table [TJ by varying the value 

of Ti. 

Fig. |4] shows the simulation results for the model R7 
series. The top and middle panels show the radial profiles of 
inflow, outflow, and net rates for the four models with dif- 
ferent ri, while the bottom panel shows the time variation 
of the luminosity of the models. The values of the outflow 
rate of the four models are almost same. Moreover, the ra- 
dius where the outflow begins to become equal to the net 
rate is also roughly same, i.e., ~ (1 — 4) x 10*r*. For the 
inflow rate at the inner boundary, models R7a and R7b are 
almost same; while model R7c is only about 10% smaller 



than model R7a and model R7d is only about 12% smaller 
than model R7c. The deviations are within the acceptable 
range. Therefore, we think model 7 series are convergent as 
we vary the radius of the inner boundaries. 

We also did the simulations with different ri for model 
R5, R6 and R8 series. The results are also listed in Ta- 
ble [l] For model R5 series, the largest discrepancy of time- 
averaged luminosity among the four models is ~ 36%. For 
model R6 and R8 series, they are 25% and 31%, respectively. 
We found that the results are similar when re-radiation is 
not included, as shown by Models K6a, K6c, and K6d. The 
change of the luminosity with Vi is complicated: the lumi- 
nosity can increase or decrease with increasing ri. Physi- 
cally, when ri increases, the attenuation of X-ray flux will 
become weaker, thus the corresponding radiation force be- 
comes stronger. But on the other hand, when X-ray radia- 
tion becomes weaker, the ionization will become weaker and 
correspondingly the line force stronger. Thus the flnal re- 
sult depends on the competition of these two effects. But 
one relatively robust result is that the outflow mainly origi- 
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Figure 3. The correlation between tiie mass outflow rate and 
Eddington ratio (F = Lacc/LEdd)- All the lines and dots except 
the three solid blue ones are taken from KP09. The three solid 
blue dots are the results when the re-radiation effect is included 
when To = 2 X lO^K. We can see that the correlation index q 
almost remains unchanged. 



nates from ~ (1 — 4) x lO^r*. Another noteworthy result is 
that the emitted luminosities shown in Fig. |4] are all super- 
Eddington, although the re-radiation effect has been taken 
into account. So the sub-Eddington puzzle is not solved for 
these simulations. 



mass flux, radial velocity, momentum flux, kinetic energy 
flux, and thermal energy flux of outflow at the outer bound- 
ary, respectively. Fig. [S] shows some additional properties of 
outflow for Model R7c. Panels (a) and (b) show the angular 
distribution of mass flux and density, respectively. We can 
see that most of the outflow occurs close to the disk surface, 
in the range of 35° < 8 < 68° . In the range oi 6 > 68° it is 
the inflowing region and there is almost no outflow. Panel 
(c) shows the radial proflle of the mass flux-weighted radial 



velocity of outflow. The typical value of Vr is 10 cm 



We 



can see that in most of the radial range, the radial veloc- 
ity does not keep increasing as we naively expect because 
of the acceleration by the radiation force. It even decreases 
somewhat at large radii. This is because of the continuous 
injection of new gas into the outflow. Panel (d) shows the 
angular distribution of Vr at ro- It monotonically increases 
from < lO^cm s~^ close to the surface of inflowing region 
at S ~ 68° to ~ 10^ cm s~^ close to the axis. The solid and 
dotted lines in Panel (e) show the radial profiles of the mo- 
mentum flux of outflow (p-w) and radiation from the central 
AGN, respectively; while the solid and dotted lines in Panel 
(f) for the radial profiles of the kinetic (Ek) and thermal en- 
ergy fluxes (-Eth) of outflow, respectively. Their deflnitions 
are: 



7r/2 



Pw{r) = 47rr / pv,. sin(^) d9 for Vr > 0, 



£fe(r) = 27vr' 



tt/2 



pv^ sin(6') d6 for Vr > 0, 



(22) 



(23) 



3.3 Dependence on gas density at the outer 
boundary 

We can see from Table[T]that the significance of re-radiative 
effects depends on the outer boundary density po- For ex- 
ample, when po — 5x lO'^^g cm~'^, the time-averaged lumi- 
nosity decreases by a factor of 2 from model K7a to model 
R7a. The outflow rate at the outer boundary increases by a 
factor of about 7 from model K7a to model R7a. But when 
Po = 10~^''g cm~^, the luminosity decreases by only ~ 35%, 
and the outflow rate increases by a factor of about 2.5 from 
model K5a to R5a. Usually we flnd that the smaller the gas 
density at the outer boundary is, the weaker the re-radiation 
effect will be. This result is easy to understand. The thick- 
ness of the accretion disk is determined by the vertical com- 
ponent of the gravitational force, gradient of gas pressure, 
and the re- radiation force. The former two forces are propor- 
tional to density while the re-radiation force is proportional 
to the square of density (refer to eauation ll4|l . Therefore, the 
disk will become thicker and more exposed to the irradiation 
from the central AGN. 



3.4 Properties of outflows 

In this section, taking models R5c, R6c, R7c and R8c as 
examples, we study the properties of the outflows. This is 
a sequence of increasing outer density and larger and lager 
rates of nominal inflow. These properties could be used to 
compare with observations. Some main properties are listed 
in Table (2] From the left to right columns, they are the 



Ethi 



4Tvr^ 



n/2 



evr sin(6') dO for Vr > 0. 



(24) 



We see that the momentum flux p™ keeps increasing out- 
wards, because more and more momentum of the radiation 
is transferred into the outflow gas. But we note that even 
at ro, the momentum flux of outflow is still several times 
smaller than that of radiation. It will be interesting to adopt 
a larger ro to see how the momentum flux of outflow can fi- 
nally become closer to that of the radiation. From Panel 
(f) we see that most of the power of outflow is in the form 
of kinetic energy rather than thermal energy. The highest 
power is reached at ro, which is ~ 1.1 x 10*^ ergs s"'^. As 
a comparison, for this model the luminosity from the cen- 
tral AGN is l.SLEdd ~ 2.3 x lO'^'^ergs s"^ This is about 
three orders of magnitude higher than the kinetic power of 
outflow. Compared to the model of KP09, the efficiency of 
transferring the radiation into kinetic power of outflow is one 
order of magnitude higher. But it is still much lower than 
~ 0.05, which is often assumed in cosmological simulations 
(see discussions in Kurosawa, Proga & Nagamine 2009). 

The value of rj^ (defined in equation I2ip is given in the 
last column of Table [l] We can see that it becomes signifi- 
cantly larger when the re- radiation effect is included, rj^ > 1. 
This implies that outflow rate is larger than the accretion 
rate. It is interesting to note that this is also the case for hot 
accretion flow (Yuan, Wu & Bu 2012) and found by Li et al. 
(2013). Of course, in that case, the mechanism of producing 
outflow is completely different. 
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Figure 4. Effects of varying tiie inner boundary r^ for model R7 series. Top and middle panels: The radial profiles of mass inflow 
rate, outflow rate, and the net rate in unit of lO^^g cm~^. Models a, b, c, and d correspond to r^ = 500r,,625r*, 1250r*, and 2500r*, 
respectively. We can see the results are roughly convergent with various r^ . Bottom panel: Time evolution of luminosity normalized by 
Eddington luminosity LEdd for model R7 series. The red, blue, green, and black lines correspond to models R7a, R7b, R7c, and R7d, 
respectively. 



4 DISCUSSION 

4.1 Effect of changing the temperature at the 
outer boundary 

So far we have fixed Tq = 2 x 10'' K. It will be interesting 
to see the effect of varying To with the adopted value To = 



2 X 10^ K, the Bondi radius is tb = GMbh/cI = 4.8 x 
10^^ cm ~ 2.2 To- If we adopt higher To but maintain the 
outer boundary unchanged, the Bondi radius will be within 
our computational domain. This is an advantage since Bondi 
radius is the place where the inflow rate is determined (Li, 
Ostriker & Sunyaev 2013). We choose three To values: 5x 10^ 
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Table 2. Properties of Outflow 



Model Mout{ro) Vr{ro) Pw{ro) E,,{ro) Eth{ro) 

(1025 g s-l) (km s-i) (10^3 g ^^ g-2) (ig40 gj.gs s-l) (lO^O ergs s"!) 



R5c 


24.28(10.67) 


800 


21.57(5.72) 


R6c 


39.10(12.13) 


1000 


43.03(16.81) 


R7c 


127.04(11.80) 


700 


98.00(5.52) 


R8c 


28.94(44.12) 


600 


23.58(64.40) 



174.16(48.80) 3.77(0.83) 

547.17(459.39) 7.56(4.09) 

1098.91(37.47) 12.31(0.63) 

1110.45(10393.40) 91.55(1206.77) 



Note: The values in brackets are the normalized standard deviations (t„ of the time series values. For the Column Vr{ro), Vr{To) is 
the time-averaged mass flux- weighted radial velocity of outflow at To (see panel c of Fig. (5)1. To avoid the influence of the outer 
boundary, we actually calculate all the quantities at ~ 6 pc. 



K, 6 X 10^ K and 10^ K. The corresponding tb of the two 
latter cases is larger than ro now. These models with new 
To are listed in the end of Table [T] (Models 26-31). 

The top panel of Fig.[6]gives the emitted luminosity Lace 
as a function of time for Models R7c (red line; To — 2 x 10® 
K), R7c-1T (blue hue; To = 5 x lO'^ K); R7c-hT (green hue; 
To = 6 X 10® K), and R7c-hT2 (black hue; To = lO'^ K). 
It is interesting to note that for model R7c-hT2 which has 
the highest Tq and smallest rs (~ 0.43ro), the light curve 
oscillates regularly. Similar behavior has been found before 
such as the feedback study in elliptical galaxies by Novak, 
Ostriker & Ciotti (2011). Such an oscillation is caused by 
radiative feedback. For model R7c-hT2, at some time the 
strong radiative heating heats the gas at a certain radius 
(j'viriai) above the local virial temperature. Then the accre- 
tion is partly suppressed since the gas becomes unbound. 
The accretion rate will thus begin to decrease, and subse- 
quently the radiation from the central AGN will decrease. 
This makes the radiative heating weaker and the temper- 
ature of the gas will decrease until below the virial value. 
So the gas supply then recovers and the accretion becomes 
active again. The timescale of such a cycle should be deter- 
mined by the accretion timescale at rviriai- From Fig. [B] we 
see that the timescale is ~ lO^^s. This is roughly the free-fall 
timescale at the outer boundary To, which is also the accre- 
tion timescale for our low angular momentum accretion. 

We would like to emphasize that for such an oscillation 
behavior to occur, it is not necessary to have rviriai ^tb- In 
fact, Yuan, Xie & Ostriker (2009) find the same oscillation 
result even when rviriai < rs when they consider the global 
Compton heating effect in hot accretion flows. The reason 
why the other three models do not have such oscillation is 
perhaps as follows. One is that the luminosity of the active 
phase from the other three models is smaller (refer to Fig. 6), 
thus the radiative heating is weaker. Another reason is that 
the specific thermal energy of the gas in Model R7c-hT2 
is the highest among the four models at the same radius, 
therefore it is the easiest to be heated to be above the virial 
temperature. 

In this context we would like to mention that the peak 
luminosity found in Yuan, Xie & Ostriker (2009) is only 
~ 2%1/Edd, much smaller than that in Model R7c-hT2, 
which is super-Eddington. One of the main reasons for the 
discrepancy is that the emitted spectrum from the innermost 
accretion flow is different. In the current work we assume a 
typical quasar spectrum, i.e., Tx ~ 8 x lO^isT and /* — 0.05. 
Such a spectrum is likely produced by a standard thin disk 
plus a corona. However, in Yuan, Xie & Ostriker (2009), 
the accretion flow is hot (such as advection-dominated ac- 



cretion flow) thus the emitted spectrum is quite different. 
The radiation temperature of the spectrum in that case is 
Tx ~ 3x l{f K and /, ~ 0.6. In this case the Compton heat- 
ing will become much stronger (e.g., Proga 2007a). In addi- 
tion, the angular distribution of radiation is also diflerent. In 
the case of hot accretion flow, the radiation is nearly spher- 
ically symmetric. This again increases the heating close to 
the equatorial plane. We have run some test simulations and 
confirmed the above eflects. And lastly, the eflectiveness of 
Compton heating is proportional to density of the flow. The 
specific angular momentum of the accretion flow in Yuan 
et al. (2009) is much larger than that in the current work. 
Therefore, for the same density, the corresponding accretion 
rate and therefore the luminosity of Yuan et al. (2009) will 
be significantly smaller than in the current work. 

4.2 Effect of X-ray optical depth 

So far we have treated the Thompson scattering opacity for 
the calculation of X-ray optical depth as if it were absorp- 
tion. However, as we have mentioned above, this is likely an 
over-estimation (we neglect X-ray absorption). Another ex- 
treme is to take the X-ray optical depth is zero. We now test 
this case. The last few models labeled with 'XO' in Table [T] 
are such models. 

Fig. [7] compares the logarithmic density and temper- 
ature contours with (model R7c, top panels) and with- 
out (model R7c-X0, bottom panels) X-ray opacity. When 
hix = 0, there will be no X-ray attenuation. In this case, 
the ionization parameter ^ increases and X-ray heating is 
enhanced especially at large radii (cf. equations I lOm If) . We 
can see from the bottom-right panel of the figure that the gas 
temperature at low latitudes increases significantly. The gas 
density near the pole decreases greatly as the X-ray pushes 
the gas away. Once the gas density becomes lower, the X- 
ray heating becomes less effective, hence the temperature 
decreases. Note in table [T] the high values of 77^, (cf. equa- 
tion I2ip in the models having tx ~ and ro > rs . The 
outflow rates exceed the net accretion rates. 

The top panel of Fig. [8] shows the light curves of var- 
ious models which have been shown in Fig. [51 except here 
Kx = 0. In Fig.[6l only the model with the highest Tq (model 
R7c-hT2) oscillates. But here with the exception of model 
R7c-1TX0 which has the lowest Tq, the light curves of all 
other three models oscillate. The discrepancy is because of 
the enhanced X-ray heating. Comparing model R7c-hT2 in 
Fig. [5] and model R7c-hT2X0 in this flgure, we can see that 
the timescale of the variability becomes smaller. This indi- 
cates that the virial radius rviriai becomes smaller. Another 
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effect is that the amplitude of oscillation of model R7c-hT2 
becomes smaller now. This is perhaps because a smaller frac- 
tion of the gas located at rviriai is heated to be unbound now. 



4.3 The luminosity observed at infinity 

The Ught curves shown in the top panels of Fig. |6] and Fig. [8] 
are the emiiterf luminosity Lace, not the observed luminosity 
at infinity (Lobs) since absorption has not been taken into 
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account. To calculate Lobs, for simplicity, we assume both 
the UV and X-ray opacities are equal to Thomson scattering 
opacity, and use our simulation data to calculate the optical 
depth along various 6 by integrating the optical depth along 
this 6 line from the inner boundary to the outer boundary. 
For the optical and UV output this simple formulation may 
lead to an overestimation of the luminosity reaching infin- 
ity. Then we can obtain the observed luminosity following 

-'^obs — ^acc6 

The bottom panels of Fig. [6] and Fig. [8] give the Lobs 
light curves. It is interesting to note that the oscillation am- 
plitude of Lobs in Fig. [6] becomes larger than Lace while in 
Fig. [8] it remains almost unchanged. This indicates that in 
the case of Fig. [6l the optical depth is larger during the in- 
active phase than that in the active phase, while in the case 
of Fig. [8] the optical depth remains almost the same in the 
two phases. Again, in both cases, the observed luminosity 
can well exceed the Eddington luminosity. 



4.4 Observational implications 

We now briefly discuss the observational implications of 
our results. The first issue is the origin of outflow ob- 



served in AGNs. They have been widely observed in the 
form of absorption lines in a variety of AGNs, includ- 
ing radio-quiet and radio loud, Seyfert and quasars, low- 
luminosity or high-luminosity ones (see Crenshaw, Kraemer 
& George 2003 for a review). Especially, recently Tombesi 
et al. (2010, 2011, 2012) detected ultrafast outflows with ve- 
locities up to (0.03 — 0.4)c via Fe K-shell absorption lines. 
They found that the ultrafast outflows are located in the in- 
terval ~ (0.0003 - 0.03) pc ~ (10^ - 10*) Ts from the black 
hole on average. The mass outflow rates are constrained be- 
tween ~ (0.01 — 0.1) M0 yr~^ and the average lower and 
upper limits on the mechanical power are Ek — 10*^'' and 
lO'**'' ergs s~^, respectively. These features are difficulty to 
explain with our simulations. In our models, the highest out- 
flow velocity is about 0.03c (see panel (d) of Fig. [5]), and the 
outflows starts from ~ 10"" r^ (see Fig. |3| , the mass outflow 
rates and outflow kinetic power at the outer boundary are 
> 4 Mq yr-i and Ek ~ (10*^ - 10*^) ergs s"^ (see Table 
[5}. However, based on these results we can not rule out the 
"line-driven" models as the origin of outflow. This is because 
our simulations concentrate only on the large scale of the ac- 
cretion flow, (0.01 ~ 7) pc, the innermost accretion flow is 
not included and substantial radiative acceleration may oc- 
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Figure 7. Comparison of the logarithmic density (left) and temperature (right) contours between Model R7c (top panels) and Model 
R7c-X0 (bottom panels; with kx = 0). The arrows represent the direction of velocity vector but not its magnitude. 



cur within n. We did adopt different inner boundary ri and 
found that the properties of outflow, including the location 
of origin, remain largely unchanged. This is likely because 
the angular momentum of the accretion flow we consider is 
rather small so there is no real disk formed in our simula- 
tion domain. On the other hand, the discrepancies between 
our simulation results and observations indicate that at least 
the ultrafast outflow must come from the innermost region 
of the accretion disk. 



5 SUMMARY 

In this paper we study the dynamics of accretion flow within 
~ 0.01 — 7 pc from the central super-massive black hole irra- 
diated by the central AGN. We focus on the outflow driven 
by the radiation via the Compton scattering and the line 
force. Compared to the previous works, the main improve- 
ment is that we include the re-radiation force, which is the 
additional force produced by the scattered photons and the 



locally produced photons. We flnd that the results change 
signiflcantly. The accretion flow now becomes thicker and is 
thus more exposed to the radiation from the central source. 
Consequently, the outflow rate measured at the outer bound- 
ary increases by about one order of magnitude. We find that 
the correlation between the Eddington ratio F and the out- 
flow rate Mout originally found in previous works still exists. 
The properties of the outflows are calculated and compared 
to observations, such as the ratio of the outflow rate and the 
accretion rate. We find that some observed outflows can not 
be explained by this model, but note that the inferred high 
ratios of Mout/ Mace occur naturally in our high Edding- 
ton ratio models. When the density of the gas at the outer 
boundary is sufficiently high, we still find super-Eddington 
accretion, even though we include the re-radiation effect and 
the outffow becomes stronger. We also investigate the effect 
of different temperatures at the outer boundary (To)- We 
find that when To is large (~ lO'^ K), the emitted luminos- 
ity from the accretion fiow oscillates significantly because of 
the radiative heating feedback. 
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One of the main aims of the present research is the sub- 
Eddington puzzle. Observations show that most AGNs are 
sub-Eddington, while on the other hand both analytical the- 
ory and numerical simulations of accretion flows have clearly 
shown that super-Eddington accretion can exist. In this pa- 
per, we hope to investigate whether the feedback at larger 
scale can solve the puzzle, especially after the re-radiation 
effect is taken into account. Unfortunately, we find that the 
super-Eddington accretion still exists (refer to Figs. [S] & HI) 
but at a reduced rate. 

There are some caveats in our study. These issues can 
be improved in the future, which may be helpful to solve the 
sub-Eddington puzzle. The first one is that we need to con- 
sider the probable large angular momentum of the accretion 
fiow. When the angular momentum is large, we must include 
viscosity and a disk will be formed. In this case, such a disk 
itself will emit strong radiation, which may help to drive an 
outfiow. Moreover, when the angular momentum is large, a 
strong outfiow will be produced intrinsically by other mech- 
anisms (Yuan, Bu & Wu 2012; Bu et al. 2013). For example, 
outflow can be produced by buoyancy or magnetocentrifugal 
force in a hydrodynamical or magnetodynamical accretion 
flow, respectively (Yuan, Bu & Wu 2012; see also Bai & 
Stone 2013). The second caveat is that in the current simu- 



lations, we only consider the radiation from within r^ and we 
do not investigate the dynamics within r^. However, strong 
outflows have been shown to exist from the innermost ac- 
cretion flow and they will interchange their momentum and 
energy with the gas. So the feedback of the outflow should 
also be included (e.g., Ostriker et al. 2010; Novak et al. 2011; 
Choi et al. 2012). All the above-mentioned factors are ex- 
pected to be able to make the outflow stronger and thus 
the inflow rate at r^ and luminosity smaller. In addition, a 
simplified approach is adopted when we deal with the in- 
teraction between the radiation, from both the central force 
and re-radiation, and the gas. Obviously a proper radia- 
tive transfer calculation is desired. And lastly, the opacity 
due to dust can also potentially greatly enhance the outfiow 
(e.g., Novak et al. 2012; Roth et al. 2012; Faucher-Giguere, 
Quataert & Hopkins 2013). 



The values oi rj^ > 1 found in our high infiow models 
are consistent with recent assessments of Moe et al. (2009) , 
Dunn et al. (2010) and the references therein that winds 
carry out more mass in bright quasars than is being accreted 
on the central black holes. 
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